TD5_Baptiste_DUFOUR

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
## Le chargement a nécessité le package : xts
## Le chargement a nécessité le package : zoo
## 
## Attachement du package : 'zoo'
## Les objets suivants sont masqués depuis 'package:base':
## 
##     as.Date, as.Date.numeric
## Le chargement a nécessité le package : TTR

Exercice 1 : Data settings

1) Data importation

fredr_set_key("6e17c533d47a0aa5de126acc5fe81ed9")

chomage_us<-fredr(
  series_id="UNRATENSA" ,
  observation_start =as.Date("1948-01-01"),
  observation_end=as.Date("2023-11-01"),
  frequency = "m" # monthly
)

2) Ploting and observing

On observe une série cyclique et saisonniaire.

Un cycle peut être comparé à une séquence d’événements qui se répètent au fil du temps. Il démarre à un point bas local de la série, atteint son point culminant, puis se termine pour débuter le cycle suivant à ce même point. Contrairement à un motif saisonnier, les cycles ne suivent pas nécessairement des intervalles de temps réguliers, et leur durée peut varier d’un cycle à l’autre.

La composante saisonnière (ou saisonnalité) est un autre motif fréquent dans les données de séries temporelles. Si elle est présente, elle représente une variation répétée dans la série, liée aux unités de fréquence de la série (par exemple, les mois de l’année pour une série mensuelle). Un exemple courant de série présentant un motif saisonnier marqué est la demande d’électricité ou de gaz naturel. Dans ces cas, le motif saisonnier est influencé par divers événements saisonniers tels que les conditions météorologiques, la saison de l’année et les heures d’ensoleillement. De plus, une série peut avoir plus d’un motif saisonnier.

Il y a une partie cyclique, on devine une saisonnalité, mais pas de tendance appparante. On ne peut pas conclure de la stationarite, la condition sur la variance n’est pas forcémement respectée. On voit que la variance dépend du temps autour de 800 observations.

On rappelle qu’une série est fiablement stationnaire si :

  • \(E(y_t)\) est indépendante de t

  • \(Var(y_t)\) est indépendante de t et aussi constante

  • \(Cov(y_t,y_s)\) dépend de (t-s) et non de t ou s.

On observe également que l’ACF décroit lentement, ce qui montre une série avec de la mémoire donc potentiellement non stationnaire. On observe des pics d’autocorrelation réguliers.

De même avec la PACF, qui a un pic très fort au lag 1 ce qui explique la décroissance lente de l’ACF, toute la mémoire dépend du premier paramètre autorégressif.

On peut donc penser avec toutes ses informations à une saisonnalité additive \(y_t=T_t+S_t+c_t+ε_t\) pour une série qui semble persistente et non stationnaire.

Nous pouvons utiliser plusieurs méthodes pour traiter la saisonnalité:

  • La régression linéaire: On génère des variables indicatrices mensuelles, des dummies et on modélise notre série comme une somme de nos paramètres de régréssion qui sont nos dummies \(y_t = t_t + s_t = β_0 + β_1z_{t1} + . . . + β_pz_{tp}\). Cette régression utilise la méthode des moindres carrés qui implique que nos observations sont indépendantes ce qui n’est pas toujours le cas.

  • Moyenne mobile: Une moyenne mobile est une série temporelle construite en prenant la moyenne de plusieurs valeurs séquentielles d’une autre série temporelle. C’est une forme de convolution mathématique. Si l’on représente la série temporelle originale par y1, …, yn, alors une moyenne mobile à deux côtés de la série temporelle est donnée par \[ z_t = \frac{1}{2q+1}\ + \sum_{j=-q}^{q} y_{t+j} \quad \text{for } t = q+1, q+2, \ldots, n-q \] De manière similaire, une moyenne mobile à un seul côté de yt est donnée par: \[ z_t = \frac{1}{q+1}\ + \sum_{j=0}^{q} y_{t-j} \quad \text{for } t = q+1, q+2, \ldots, n \] Les moyennes mobiles sont utilisées de deux manières principales :

  1. Les moyennes mobiles à deux côtés (pondérées) sont utilisées pour “lisser” une série temporelle afin d’estimer ou mettre en évidence la tendance sous-jacente ;
  2. Les moyennes mobiles à un seul côté (pondérées) sont utilisées comme méthodes simples de prévision pour les séries temporelles.
  • La différenciation: Les tendances peuvent être éliminées en différenciant les données. Considérons des données avec une tendance linéaire modélisée par Xt = β0 + β1t + et où et est un processus stationnaire. La différence saisonnière consiste à prendre la différence entre une observation et l’observation précédente de la même saison (mois, trimestre, etc.).

\[y_t' = y_t - y_{t-m}\]

Ici m est le nombre de saisons. Ces différences sont également appelées “différences de retard” car nous soustrayons l’observation après un retard de m périodes. Si les données différenciées saisonnièrement semblent être un bruit blanc, alors un modèle approprié pour les données d’origine est :

\[y_t' - y_{t-m} = ε_t\]

Ici nous allons essayer d’enlever la saisonnalité en utilisant une moyenne mobile ou avec la fonction decompose qui va nous permettre de décomposer notre série en une composante saisonnale et une tendance afin de pouvoir la soustraire à notre série de base.

3) Application du filtre

Test de la moyenne mobile

Test de la fonction decompose

On peut observer que la moyenne mobile semble avoir mieux réduit la saisonnalité mais ne semble pas non plus l’avoir éliminé.

4) Test du filtre

L’ACF décroit toujours lentement et le premier coef de la PACF n’a pas changé. On peut donc penser que la série n’est toujours pas stationnaire. On va vérifier cela avec un test de stationnarité, ici le test ADF.

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression none 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.90260 -0.02875  0.00272  0.03239  0.90959 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## z.lag.1    -0.0004364  0.0003917  -1.114    0.266    
## z.diff.lag  0.8435490  0.0179287  47.050   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0699 on 897 degrees of freedom
## Multiple R-squared:  0.7117, Adjusted R-squared:  0.711 
## F-statistic:  1107 on 2 and 897 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -1.1141 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau1 -2.58 -1.95 -1.62

La statistique du test est de -1.1141, ce qui n’est pas inférieur à la valeur critique pour 5% de -1.95, on ne peut donc pas rejeter le test et donc on conserve le fait que la série n’est pas stationnaire.

On ne peut donc pas modéliser cette série sous la forme d’un ARMA(p,q).

Exercice 2 : Unit Root tests

On peut constater que les données ajustées saisonnièrement (ou filtrées) montrent toujours un schéma cyclique. Nous devons caractériser ce schéma, c’est-à-dire déterminer si la série temporelle est stationnaire ou non, et quelle est la nature de la non-stationnarité le cas échéant.

1) Explication du test ADF

Le test Augmented Dickey-Fuller est utile pour déterminer la stationnarité d’une série. On pose π = ϕˆ− 1

  • Hyptohèse du test: Ho: π=0, la série est non stationnaire

H1: π<0, la série est stationnaire.

Le test ADF est basé sur les trois équations:

\[▽y_t= πy_{t-1}+ \sum_{i=2}^{p} \gamma_i \Delta y_{t-i+1} + u_t \]

\[▽y_t= β_1 + πy_{t-1}+ \sum_{i=2}^{p} \gamma_i \Delta y_{t-i+1} + u_t \]

\[▽y_t= β_1 + β_2t + πy_{t-1}+ \sum_{i=2}^{p} \gamma_i \Delta y_{t-i+1} + u_t \]

La différence entre ces équations de test concerne l’inclusion des coefficients β1 et β2, où, sous l’hypothèse nulle, la première équation se réfère à un modèle de marche aléatoire pure, la deuxième équation inclut un terme d’intercept ou de drift, et la dernière équation inclut à la fois un drift et une tendance linéaire dans le temps. Dans chaque cas, le paramètre d’intérêt est π, et si π = 0, alors le processus yt contient une racine unitaire.

La comparaison du t-statistique calculé avec les valeurs critiques des tables de Dickey-Fuller détermine si nous devons rejeter ou non l’hypothèse nulle, H0 : π = 0

Nous allons maintenant montrer un exemple de ce test avec une marche aléatoire.

yt=cumsum(rnorm(100))
test=ur.df(yt,type="none",selectlags="BIC") 
summary(test)
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression none 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.3194 -0.4712  0.1055  0.6930  2.1743 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## z.lag.1     0.0005087  0.0295554   0.017    0.986
## z.diff.lag -0.0160854  0.1078435  -0.149    0.882
## 
## Residual standard error: 0.9294 on 96 degrees of freedom
## Multiple R-squared:  0.000241,   Adjusted R-squared:  -0.02059 
## F-statistic: 0.01157 on 2 and 96 DF,  p-value: 0.9885
## 
## 
## Value of test-statistic is: 0.0172 
## 
## Critical values for test statistics: 
##      1pct  5pct 10pct
## tau1 -2.6 -1.95 -1.61

On observe une valeur de la statistique bien au dessus de la valeur critique ce qui est logique comme cela représente une marche aléatoire donc une série non stationnaire.

Maintenant, testons avec la différentiation de cette série.

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression none 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.12971 -0.56909  0.02041  0.72125  2.30292 
## 
## Coefficients:
##            Estimate Std. Error t value Pr(>|t|)    
## z.lag.1    -1.10343    0.14460  -7.631 1.79e-11 ***
## z.diff.lag  0.09061    0.10132   0.894    0.373    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9167 on 95 degrees of freedom
## Multiple R-squared:  0.5139, Adjusted R-squared:  0.5037 
## F-statistic: 50.22 on 2 and 95 DF,  p-value: 1.316e-15
## 
## 
## Value of test-statistic is: -7.6307 
## 
## Critical values for test statistics: 
##      1pct  5pct 10pct
## tau1 -2.6 -1.95 -1.61

On observe ici une statistique bien inférieure à celle de la valeur critique donc cela montre bien une série stationnaire.

2) Différentiation de la série

Nous avons précedemment réaliser le test ADF sur notre série, ce qui a montré la non stationnarité. Nous allons maintenant essayer de différencier notre série pour la rendre stationnaire.

chomage_diff=diff(chomage_ma,lag=12,differences=1)
ts.plot(chomage_diff, main = "Chomage")

test=ur.df(chomage_diff,type="none",selectlags="BIC") 
summary(test)
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression none 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.90782 -0.02477 -0.00320  0.01989  0.93916 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## z.lag.1    -0.024552   0.002027  -12.11   <2e-16 ***
## z.diff.lag  0.922563   0.012777   72.20   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.07072 on 885 degrees of freedom
## Multiple R-squared:  0.8561, Adjusted R-squared:  0.8558 
## F-statistic:  2633 on 2 and 885 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -12.1143 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau1 -2.58 -1.95 -1.62

On observe ici une statistique de test de -12 largement inférieur à la valeur critique donc on rejette Ho et on considère la série comme stationnaire.

3) Degré d’intégration de la série

Le degré d’intégration d’une série est le nombre de différences nécessaire pour rendre la série stationnaire.

On observe que après avoir différencier la série à l’ordre 1, la série est stationnaire.

On en déduit donc que le degré d’intégration de la série est de 1.

4) Test de Philips and Perron

Le test de Philips et Perron est semblable au test de Augmented Dickey Fuller par ses hypothèses mais est un test plus paramétrique et un correctif est appliquée sur la variance estimée des résidus.

pp_test <- ur.pp(chomage_diff)
summary(pp_test)
## 
## ################################## 
## # Phillips-Perron Unit Root Test # 
## ################################## 
## 
## Test regression with intercept 
## 
## 
## Call:
## lm(formula = y ~ y.l1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.59644 -0.04876  0.00757  0.08203  1.02739 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.002599   0.006228  -0.417    0.677    
## y.l1         0.985615   0.005297 186.061   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1856 on 886 degrees of freedom
## Multiple R-squared:  0.975,  Adjusted R-squared:  0.975 
## F-statistic: 3.462e+04 on 1 and 886 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic, type: Z-alpha  is: -58.5159 
## 
##          aux. Z statistics
## Z-tau-mu           -0.2372

On observe que la p-value est très faible, ce qui indique que l’on rejette Ho donc la série est stationnaire. Cela concorde avec le test de ADF réalisé précedemment.

5) Test de KPSS

Le test de KPSS est également utilisé pour déterminer si une série est stationnaire mais varie par ses hypothèses.

  • Ho: La série est stationnaire
  • H1: La série n’est pas stationnaire.

En comparant la valeur du test-statistic avec les valeurs critiques :

Si la valeur du test-statistic est inférieure à la valeur critique correspondant à un certain niveau de signification, on ne peut pas rejeter l’hypothèse nulle. Cela suggère que la série temporelle est stationnaire.

Si la valeur du test-statistic est supérieure à la valeur critique, on peut rejeter l’hypothèse nulle, indiquant que la série temporelle n’est pas stationnaire.

kpss_test <- ur.kpss(chomage_ma,type='tau')
summary(kpss_test)
## 
## ####################### 
## # KPSS Unit Root Test # 
## ####################### 
## 
## Test is of type: tau with 6 lags. 
## 
## Value of test-statistic is: 0.6353 
## 
## Critical value for a significance level of: 
##                 10pct  5pct 2.5pct  1pct
## critical values 0.119 0.146  0.176 0.216

Pour notre série trouvé dans l’exercice 1 avant différentiation, la valeur du test-statistic est supérieure aux valeurs critiques. On peut donc rejeter l’hypothèse nulle de stationnarité. Cela suggère que la série temporelle peut être considérée comme non stationnaire.

6) Degré d’intégration KPSS

Nous allons maintenant procéder au même test mais après différentiation à l’ordre 1.

kpss_test <- ur.kpss(chomage_diff,type='tau')
summary(kpss_test)
## 
## ####################### 
## # KPSS Unit Root Test # 
## ####################### 
## 
## Test is of type: tau with 6 lags. 
## 
## Value of test-statistic is: 0.0448 
## 
## Critical value for a significance level of: 
##                 10pct  5pct 2.5pct  1pct
## critical values 0.119 0.146  0.176 0.216

Ici la valeur du test-statistic est inférieure aux valeurs critiques. On ne peut donc pas rejeter l’hypothèse nulle de stationnarité. Cela suggère que la série temporelle peut être considérée comme stationnaire.

Cela valide notre hypothèse que le degré d’intégration de la série est 1.

Exercice 3 : Modeling

1) Proposer un model ARMA(p,q)

Il existe deux types de méthodes pour identifier l’ordre du processus ARMA. Le premier repose sur les informations provenant de l’ACF et du PACF (autocorrélation partielle). À titre de rappel, l’ACF est utilisée pour évaluer l’ordre d’un processus MA(q), tandis que la PACF est employée pour déterminer l’ordre d’un processus AR(p). La deuxième méthode consiste à utiliser des critères d’information.

Un critère d’information mesure la qualité d’ajustement du modèle en attribuant un coût informationnel au nombre de paramètres à estimer.

Les critères d’information intègrent deux facteurs :

  1. Un terme qui est une fonction de la somme résiduelle des carrés ;
  2. Un terme qui pénalise la perte de degrés de liberté due à l’ajout de paramètres supplémentaires.

L’objectif est donc de choisir le nombre de paramètres qui minimise la valeur du critère d’information.

Le critère d’information d’Akaike (AIC), le critère d’information bayésien de Schwarz (SBIC) et le critère d’information de Hannan-Quinn (HQIC) sont définis par les relations suivantes :

\[AIC=2k-2ln(\hat L)\]

\[BIC=k*ln(n)-2ln(\hat L)\]

\[HQIC=-2\hat L{max}+2k*ln(ln(n))\]

et où \(\hat L\) represente la log vraissemblance, k le nombre de paramètres et n le nombre d’observations.

Méthode empirique : Nous procédons par balayage, en définissant une plage de valeurs pour p et q, et en calculant les critères d’information pour toutes les combinaisons possibles. Nous sélectionnons ensuite la paire (p, q) qui minimise les critères d’information.

Afin de définir cette plage, nous allons regarder l’ACF et la PACF.

On observe que l’ACF est en dessous du seuil pour un lag d’environ 4 et la PACF pour un lag d’environ 1 ou 2.

On va donc tester les critères d’information pour des valeurs de p et q allant de 1 à 4.

##           [,1]      [,2]      [,3]      [,4]
## [1,] -2165.561 -2164.095 -2169.726 -2167.763
## [2,] -2164.779 -2174.306 -2168.035 -2177.046
## [3,] -2170.036 -2241.463 -2238.827 -2203.639
## [4,] -2170.499 -2208.377 -2243.324 -2251.075
##           [,1]      [,2]      [,3]      [,4]
## [1,] -2151.191 -2144.935 -2145.776 -2139.022
## [2,] -2145.619 -2150.355 -2139.294 -2143.515
## [3,] -2146.085 -2212.722 -2205.296 -2165.318
## [4,] -2141.758 -2174.846 -2205.003 -2207.964
##           [,1]      [,2]      [,3]      [,4]
## [1,] -2160.068 -2156.771 -2160.571 -2156.777
## [2,] -2157.455 -2165.151 -2157.049 -2164.229
## [3,] -2160.881 -2230.477 -2226.010 -2188.991
## [4,] -2159.513 -2195.561 -2228.677 -2234.596
## [1] "Optimal AIC - p: 4 q: 4"
## [1] "Optimal BIC - p: 3 q: 2"
## [1] "Optimal QIC - p: 4 q: 4"

On trouve des critères d’information optimals pour p valant 3 et q valant 2. On va donc estimer un model ARMA(3,2)

## 
## Call:
## arima(x = chomage_ts, order = c(3, 1, 2))
## 
## Coefficients:
##          ar1      ar2      ar3      ma1     ma2
##       0.9528  -0.9534  -0.0465  -0.9927  0.9975
## s.e.  0.0332   0.0332   0.0332   0.0043  0.0103
## 
## sigma^2 estimated as 0.2701:  log likelihood = -699.24,  aic = 1410.48

2) Tests de qualité des résidus

Pour assurer la validité du model, on va devoir réaliser plusieurs tests de qualité.

Distribution Normale des Résidus :

L’affirmation selon laquelle les résidus suivent une distribution normale est cruciale pour la validité du modèle de régression linéaire. Cette hypothèse garantit que les erreurs sont impartiales et constantes à travers différents niveaux de la variable prédictive. L’histogramme offre une représentation visuelle de la distribution des résidus. Un histogramme en forme de cloche, ressemblant à une distribution normale, suggère que les erreurs du modèle sont réparties symétriquement autour de zéro, renforçant la fiabilité du modèle. Le graphique quantile-quantile (Q-Q) est un autre outil de diagnostic qui compare la distribution des résidus à une distribution normale attendue. Un motif diagonal dans le graphique Q-Q indique que les résidus s’alignent étroitement avec une distribution normale. Le test de Shapiro est également utile pour déterminer si les résidus suivent une distribution normale.

Absence d’Autocorrélation :

L’autocorrélation dans les résidus implique que les erreurs à un moment donné sont corrélées avec les erreurs à d’autres moments, indiquant potentiellement que le modèle n’a pas capturé de manière adéquate les motifs temporels. Le test de Ljung-Box évalue l’absence d’autocorrélation dans les résidus sur différents intervalles de décalage. Une valeur p non significative dans ce test indique que les autocorrélations ne sont pas statistiquement différentes de zéro, soutenant l’absence d’autocorrélation. Le graphique de la fonction d’autocorrélation (ACF) offre une représentation visuelle des autocorrélations à différents intervalles de décalage. Une diminution rapide des valeurs d’autocorrélation suggère que les résidus n’expriment pas de motif significatif au fil du temps.

Nous connaissons différents tests d’autocorrélation:

Test de Ljung-Box : La fonction d’autocorrélation (ACF) et la fonction d’autocorrélation partielle (PACF) sont des outils qualitatifs utiles pour évaluer la présence d’autocorrélation à des retards individuels. Le test Q de Ljung-Box est une approche plus quantitative pour tester l’autocorrélation à plusieurs retards simultanément. L’hypothèse nulle de ce test est que les premières m autocorrélations sont conjointement nulles, H0 : ρ1 = ρ2 = . . . , = ρm = ρm = 0 En d’autres termes :

H0 : Il n’y a pas d’autocorrélation dans les données H1 : Il existe une autocorrélation significative

La formule de ce test est :

\[ Q = n(n+2) \sum_{k=1}^{h} \frac{\hat{\rho}_k^2}{n-k} \]

Si la p-value > 5%, nous ne pouvons pas rejeter H0 à 95%.

Si la p-value < 5%, nous pouvons rejeter H0 à 95% et affirmer qu’il existe une autocorrélation significative dans les données.

Le test de Box-Pierce : test statistique utilisé pour évaluer la présence d’autocorrélation dans une série temporelle. Il s’agit d’un cas spécifique du test de Ljung-Box, plus général. Les deux tests sont couramment utilisés pour le diagnostic de modèles dans l’analyse des séries temporelles. La statistique de test de Box-Pierce est calculée comme :

\[ Q = n \sum_{k=1}^{h} \hat{\rho}_k^2 \]

Les hypothèses pour le test de Box-Pierce sont identiques à celles du test de Ljung-Box :

H0 : Il n’y a pas d’autocorrélation dans les données H1 : Il existe une autocorrélation significative

Le test de Durbin-Watson : Alternative au test de Ljung-Box. Il s’agit d’une statistique de test utilisée pour détecter la présence d’autocorrélation au décalage 1 dans les résidus (erreurs de prédiction) d’une analyse de régression. Durbin et Watson ont proposé un test basé sur une mesure de la distance entre les erreurs en i et en i − 1.

Si \(e_t\) Les résidus donnés par \(e_t=ρe_{t-1}+vt\), le test de Durbin-Watson est:

\[ d = \frac{\sum_{t=2}^{n} (ê_t - ê_{t-1})^2}{\sum_{t=1}^{n} ê_t^2} \] où T est le nombre d’observations.

Cette statistique peut être exprimée approximativement pour de grandes valeurs de T en fonction du coefficient d’autocorrélation des résidus, noté généralement ρ. \(d ∼ 2(1 − ρ)\)

Le test peut donc être interprété comme suit :

Si ρ = 0 (c’est-à-dire aucune autocorrélation), alors d est proche de 2. Si ρ = 1 (c’est-à-dire autocorrélation positive), alors d est proche de 0. Si ρ = −1 (c’est-à-dire autocorrélation négative), alors d est proche de 4.

Ainsi, la statistique de Durbin-Watson varie de 0 à 4. Une valeur autour de 2 suggère l’absence d’autocorrélation, tandis que des valeurs significativement inférieures à 2 indiquent une autocorrélation positive, et des valeurs significativement supérieures à 2 indiquent une autocorrélation négative.

Homoscédasticité

Afin de se rendre compte si les résidus sont homoscédastiques, à variance constante, on va prendre les résidus et les élevés au ². On aura donc une approximation de la volatilité. Ensuite on regarde l’ACF ou la PACF. Une deuxième méthode serait d’estimer un model AR sur sur ε² et regarder les paramètres du modèle. Si les résidus sont homoscedastiques, alors les paramètres devraient être environ de 0, si les résidus sont heteroscedastiques, les paramètres seraient différents de 0.

Stationnarité

Nous devons également regarder la stationnarité des résidus. Pour cela, nous allons réaliser le test d’ADF, présenté précedemment.

Si toutes ces conditions sont respectées, on peut ainsi valider notre model.

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals_model
## W = 0.66984, p-value < 2.2e-16

On observe que notre histogramme ne semble pas correspondre exactement à celui d’une loi normale mais les quantiles de nos résidus semblent s’approcher de ceux d’une loi normale. De même, le test de shapiro affiche une p-value très faible donc on rejette le fait que les résidus se comportent comme une loi normale.

## 
##  Box-Ljung test
## 
## data:  residuals_model
## X-squared = 3.3739, df = 5, p-value = 0.6426

On observe une p-value supérieure à 5% pour le test de Ljung-Box, on peut donc conserver l’hypothèse d’absence d’autocorrélation de nos résidus, cela est également vérifier avec l’ACF qui se rapproche très rapidement de 0.

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression none 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.0018 -0.2530 -0.0358  0.1981 10.1637 
## 
## Coefficients:
##            Estimate Std. Error t value Pr(>|t|)    
## z.lag.1    -1.04724    0.04693 -22.317   <2e-16 ***
## z.diff.lag  0.04489    0.03315   1.354    0.176    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5196 on 907 degrees of freedom
## Multiple R-squared:  0.5025, Adjusted R-squared:  0.5014 
## F-statistic: 458.1 on 2 and 907 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -22.3167 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau1 -2.58 -1.95 -1.62

on observe que la statistique du test de ADF est très faible, environ -22, largement inférieure à la valeur critique, on peut donc rejeter Ho et dire que notre série des résidus est stationnaire.

## 
## Call:
## ar(x = residuals_2, order.max = 1)
## 
## 
## Order selected 0  sigma^2 estimated as  11.8

En estimant un model AR(1) sur les résidus au ², on se rend compte que le coefficient est très proche de 0, on peut donc considérer les résidus comme homoscédastiques.

Malgré le fait que nos résidus ne se comportent pas exactement comme une loi normale, il semble y avoir une absence d’autocorrélation, de la stationnarité, et de l’homoscédasticité, on peut donc valider notre model.

Exercice 4 : Estimating an ARIMA(p,d,q)

Comme observé pendant la séance, les modèles ARIMA ont été conçus pour traiter les séries temporelles non stationnaires. Nous proposons d’utiliser ce type de spécification pour modéliser le cours de l’action Johnson & Johnson depuis 1997 jusqu’à présent, sur une base mensuelle.

On commence par charger les données:

## [1] "JNJ"

On observe une série qui comporte plusieurs tendances avec notamment une potentielle coupure vers 200 observations mais pas vraiment de saisonnalité.

1) Degré d’intégration

Afin de déterminer le degré d’intération de la série des prix de Johnson & Johnson, nous allons réaliser différents tests d’ADF pour des intégrations différentes. Lorsque la série sera considéré comme stationnaire à un degré I, nous pourrons alors considéré que la série est intégrable à l’ordre I.

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression trend 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -15.613  -1.514   0.133   1.593  16.202 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -0.07810    0.42447  -0.184   0.8541  
## z.lag.1     -0.02308    0.01201  -1.922   0.0555 .
## tt           0.01279    0.00589   2.171   0.0306 *
## z.diff.lag  -0.04426    0.05597  -0.791   0.4296  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.701 on 318 degrees of freedom
## Multiple R-squared:  0.01745,    Adjusted R-squared:  0.008183 
## F-statistic: 1.883 on 3 and 318 DF,  p-value: 0.1323
## 
## 
## Value of test-statistic is: -1.9223 3.2463 2.385 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau3 -3.98 -3.42 -3.13
## phi2  6.15  4.71  4.05
## phi3  8.34  6.30  5.36

On observe une valeur de la statistique de -1.92 ce qui est supérieur à -3.42 la valeur critique de tau3, on ne peut donc pas rejeter Ho et on considère la série comme non stationnaire.

Nous allons maintenant réaliser de nouveau ce test mais en différenciant la série à l’ordre 1.

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression trend 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.5792  -1.3218   0.0534   1.5058  13.3087 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.142005   0.406668   0.349    0.727    
## z.lag.1     -1.312146   0.079475 -16.510  < 2e-16 ***
## tt           0.002676   0.002184   1.225    0.221    
## z.diff.lag   0.243045   0.054791   4.436 1.27e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.617 on 317 degrees of freedom
## Multiple R-squared:  0.5552, Adjusted R-squared:  0.551 
## F-statistic: 131.9 on 3 and 317 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -16.5102 90.8677 136.2967 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau3 -3.98 -3.42 -3.13
## phi2  6.15  4.71  4.05
## phi3  8.34  6.30  5.36

On observe maintenant une statistique de test de -16, ce qui est largement inférieur aux valeurs critiques. On peut donc rejeter Ho et on considère la série stationnaire.

Après une différenciation à l’ordre 1, la série est devenu stationnaire, le degré d’intégration de la série est donc de 1.

2) Déterminer un model ARIMA

Nous allons maintenant estimer les coefficients p et q grâce aux critères d’informations.

## Warning in log(s2): Production de NaN

## Warning in log(s2): Production de NaN
##          [,1]     [,2]     [,3]     [,4]     [,5]     [,6]     [,7]
## [1,] 2003.585 1771.922 1772.418 1750.647 1751.297 1753.128 1751.406
## [2,] 1944.719 1773.148 1766.137 1751.319 1748.706 1745.803 1743.501
## [3,] 1891.709 1756.655 1751.118 1748.462 1753.587 1743.981 1744.420
## [4,] 1857.196 1754.325 1753.044 1747.231 1748.868 1741.887 1742.230
## [5,] 1833.807 1754.985 1747.692 1748.804 1749.638 1741.448 1743.414
## [6,] 1833.884 1756.639 1749.665 1750.126 1742.334 1743.392 1745.439
## [7,] 1785.448 1742.426 1740.610 1739.406 1738.032 1739.522 1732.963
##          [,1]     [,2]     [,3]     [,4]     [,5]     [,6]     [,7]
## [1,] 2014.918 1787.033 1791.307 1773.313 1777.740 1783.349 1785.405
## [2,] 1959.829 1792.036 1788.803 1777.763 1778.927 1779.802 1781.278
## [3,] 1910.598 1779.321 1777.562 1778.683 1787.586 1781.757 1785.975
## [4,] 1879.862 1780.769 1783.266 1781.229 1786.644 1783.442 1787.562
## [5,] 1860.251 1785.207 1781.691 1786.581 1791.192 1786.780 1792.523
## [6,] 1864.105 1790.638 1787.442 1791.680 1787.666 1792.502 1798.326
## [7,] 1819.447 1780.203 1782.164 1784.738 1787.142 1792.409 1789.628
##          [,1]     [,2]     [,3]     [,4]     [,5]     [,6]     [,7]
## [1,] 2008.109 1777.954 1779.958 1759.695 1761.852 1765.192 1764.978
## [2,] 1950.751 1780.688 1775.185 1761.875 1760.770 1759.375 1758.581
## [3,] 1899.249 1765.703 1761.674 1760.526 1767.159 1759.061 1761.008
## [4,] 1866.244 1764.881 1765.108 1760.803 1763.947 1758.475 1760.326
## [5,] 1844.363 1767.049 1761.264 1763.884 1766.226 1759.544 1763.018
## [6,] 1845.948 1770.211 1764.745 1766.714 1760.430 1762.996 1766.551
## [7,] 1799.020 1757.506 1757.197 1757.502 1757.636 1760.634 1755.583
## [1] "Optimal AIC - p: 6 q: 6"
## [1] "Optimal BIC - p: 0 q: 3"
## [1] "Optimal QIC - p: 6 q: 6"

On observe que BIC est optimal pour p=0 et q=3. On essaye de vérifier ces valeurs avec l’ACF et la PACF.

On observe que la fonction d’autocorrelation a des valeurs différentes de 0 jusqu’au lag 3 environ, ce qui peut se rassembler à un q=3.

La fonction d’autocorrelation partielle elle semble aller vers 0 vers le lag 1.

Les observations de l’ACF et la PACF sont globalement les mêmes que celles avec les critères d’information. On va donc générer un model ARMA(0,3) pour modéliser notre série de prix.

3) Création du model et affichage

predicted_model <- arima(jnj_ts, order = c(0, 1, 3))
predicted_model
## 
## Call:
## arima(x = jnj_ts, order = c(0, 1, 3))
## 
## Coefficients:
##           ma1      ma2      ma3
##       -0.0673  -0.2239  -0.0275
## s.e.   0.0555   0.0591   0.0641
## 
## sigma^2 estimated as 13.13:  log likelihood = -874.27,  aic = 1756.54
jnj_arma=fitted(predicted_model)


plot(jnj_ts, type = "l", col = "blue", ylab = "Stock Price", xlab = "Time",main="Observed Vs Predicted Values")
lines(fitted(predicted_model), col = "red")
legend("topleft", legend = c("Observed", "Predicted Values"), col = c("blue", "red", "green"), lty = 1)

On observe que le model semble suivre notre courbe des prix. Le model ARMA(0,3) semble donc potentiellement être bien adapté pour notre série.

4) Tests de qualité des résidus

Comme dans l’exercice précédent, nous allons procéder aux tests de qualité: distribution normale des résidus, absence d’autocorrélation, homoscédasticité et stationnarité.

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals_model
## W = 0.92495, p-value = 1.104e-11

On observe que notre histogramme semble correspondre à celui d’une loi normale mais les quantiles extrèmes de nos résidus s’éloignent de ceux d’une loi normale. De même, le test de shapiro affiche une p-value très faible donc on rejette le fait que les résidus se comportent comme une loi normale.

## 
##  Box-Ljung test
## 
## data:  residuals_model
## X-squared = 2.8453, df = 5, p-value = 0.7238

On observe une p-value supérieure à 5% pour le test de Ljung-Box, on peut donc conserver l’hypothèse d’absence d’autocorrélation de nos résidus, cela est également vérifier avec l’ACF qui se rapproche très rapidement de 0.

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression none 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.2367  -0.8548   0.6424   2.0735  14.7871 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## z.lag.1    -1.007018   0.079250 -12.707   <2e-16 ***
## z.diff.lag  0.006599   0.056115   0.118    0.906    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.641 on 320 degrees of freedom
## Multiple R-squared:  0.5002, Adjusted R-squared:  0.4971 
## F-statistic: 160.1 on 2 and 320 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -12.7069 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau1 -2.58 -1.95 -1.62

on observe que la statistique du test de ADF est très faible, environ -12, largement inférieure à la valeur critique, on peut donc rejeter Ho et dire que notre série des résidus est stationnaire.

## 
## Call:
## ar(x = residuals_2, order.max = 1)
## 
## Coefficients:
##      1  
## 0.1657  
## 
## Order selected 1  sigma^2 estimated as  862.3

En estimant un model AR(1) sur les résidus au ², on se rend compte que le coefficient est très proche de 0, on peut donc considérer les résidus comme homoscédastiques.

Malgré le fait que nos résidus ne se comportent pas exactement comme une loi normale, il semble y avoir une absence d’autocorrélation, de la stationnarité et de l’homoscédasticité, on peut donc valider notre model.

4) Prédiction

## Warning in forecast.Arima(predicted_model, h = 3, xlim = ): The non-existent
## xlim arguments will be ignored.

##               80%      95%
## Jan 2024 151.0201 148.5614
## Feb 2024 149.0508 145.6888
## Mar 2024 148.2352 144.4485
##               80%      95%
## Jan 2024 160.3089 162.7675
## Feb 2024 161.7527 165.1147
## Mar 2024 162.5420 166.3288

Les tableaux affichent les valeurs inférieures de l’interval de confiance puis les valeurs supérieures.

Exercice 5 : Unit root test another one

1) Présentation du test

Dans cet exercice, nous proposons d’introduire un nouveau test de racine unitaire : le test de Zivot and Andrews (1992). Une des faiblesses des tests de racine unitaire ADF est leur potentiel à confondre les ruptures structurelles dans la série comme des preuves de non-stationnarité. En d’autres termes, ils peuvent échouer à rejeter l’hypothèse de racine unitaire si la série présente une rupture structurelle.

Le test de rupture structurelle endogène de Zivot and Andrews (1992) est un test séquentiel qui utilise l’échantillon complet et utilise une variable factice différente pour chaque date possible de rupture. La date de rupture est sélectionnée là où la statistique t du test ADF de la racine unitaire est minimale (la plus négative). Par conséquent, une date de rupture sera choisie là où les preuves sont les moins favorables pour l’hypothèse de racine unitaire.

Les tests de Zivot-Andrews (1992) énoncent l’hypothèse nulle selon laquelle la série a une racine unitaire avec des ruptures structurelles contre l’hypothèse alternative selon laquelle elles sont stationnaires avec des ruptures. Nous rejetons l’hypothèse nulle si la statistique t est inférieure à la valeur critique tabulée.

2) Génération de 3 marches aléatoires

Nous allons générer trois marches aléatoires. La première étant une pure marche aléatoire, la deuxième une marche aléatoire avec une rupture au milieu et la troisième, une marche aléatoire avec une rupture au milieu et dans la tendance. Une quatrième série accentuant les ruptures sera générer.

e=rnorm(1000)
trd=1:1000
s=c(rep(0,499),rep(1,501))
s2=c(rep(0,299),rep(1,701))

y1=cumsum(e)
y2=s*30+cumsum(e)
y3=0.1*trd+s2*25+cumsum(e)
y4=trd*s*0.1+s2*50+cumsum(e)

par(mfrow=c(2, 2))
plot(y1, type="l", main="Marche aléatoire pure", col="black")
plot(y2, type="l", main="Marche aléatoire avec rupture au milieu", col="red")
plot(y3, type="l", main="Marche aléatoire avec rupture au milieu et dans la tendance", col="blue")
plot(y4, type="l", col="yellow")

Afin de réaliser ces marches alatoires, nous avons créer un vecteur de bruits blancs gaussiens ainsi que un vecteur représentant une tendance trd et deux vecteurs qui vont créer des ruptures, s crée une rupture vers 499 et s2 vers 299.

On observe bien sur le graphique une rupture pour y2 en 500 et une pour y3 vers 299 et une dans la trend. On a créer également une série y4 qui accentue cette rupture dans la trend et avant afin de la comparer à y3.

3) Implémentation du test de Zivot et Andrews

Nous allons maintenant réaliser le test de Zivot et Andrews pour nos différentes séries. Pour cela, nous avons créer une fonction qui va nous permettre de collecter les informations importantes dans une table afin de pouvoir les visualiser. Nous récupérons ainsi la valeur de la statistique, les valeurs critiques et le point de rupture potentiel.

zivot_andrews_test <- function(serie, serie_name,model) {
  za.test <- ur.za(serie, model = model)
  print(summary(za.test))
  plot(za.test)
  test_statistic <- za.test@teststat
  critical_values <- za.test@cval
  rupture=za.test@bpoint
  
  # Create a summary table
  summary_table <- data.frame(
    Series = serie_name,
    `Test Statistique` = test_statistic,
    `Valeur critique 1%` = critical_values[1],
    `Valeur critique 5%` = critical_values[2],
    `Valeur critique 10%` = critical_values[3],
    `Point de rupture`=rupture
  )
  
  return(summary_table)
}

# Perform the test for each series
summary_y1 <- zivot_andrews_test(y1, "y1","both")
## 
## ################################ 
## # Zivot-Andrews Unit Root Test # 
## ################################ 
## 
## 
## Call:
## lm(formula = testmat)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.91878 -0.66127  0.00917  0.63897  3.12541 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.4191857  0.1194352  -3.510 0.000469 ***
## y.l1         0.9652107  0.0072061 133.943  < 2e-16 ***
## trend        0.0014650  0.0004049   3.618 0.000312 ***
## du          -0.7281646  0.1770289  -4.113 4.22e-05 ***
## dt           0.0008565  0.0004568   1.875 0.061062 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9943 on 994 degrees of freedom
##   (1 observation effacée parce que manquante)
## Multiple R-squared:  0.9866, Adjusted R-squared:  0.9866 
## F-statistic: 1.831e+04 on 4 and 994 DF,  p-value: < 2.2e-16
## 
## 
## Teststatistic: -4.8277 
## Critical values: 0.01= -5.57 0.05= -5.08 0.1= -4.82 
## 
## Potential break point at position: 503

summary_y2 <-zivot_andrews_test(y2, "y2","both")
## 
## ################################ 
## # Zivot-Andrews Unit Root Test # 
## ################################ 
## 
## 
## Call:
## lm(formula = testmat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.0363 -0.7074  0.0100  0.5939 30.5796 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.6893679  0.1665301  -4.140 3.77e-05 ***
## y.l1         0.9397564  0.0098129  95.767  < 2e-16 ***
## trend        0.0023379  0.0005695   4.105 4.37e-05 ***
## du           0.9299570  0.2152830   4.320 1.72e-05 ***
## dt           0.0005729  0.0006383   0.898     0.37    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.397 on 994 degrees of freedom
##   (1 observation effacée parce que manquante)
## Multiple R-squared:  0.995,  Adjusted R-squared:  0.9949 
## F-statistic: 4.914e+04 on 4 and 994 DF,  p-value: < 2.2e-16
## 
## 
## Teststatistic: -6.1392 
## Critical values: 0.01= -5.57 0.05= -5.08 0.1= -4.82 
## 
## Potential break point at position: 498

summary_y3 <- zivot_andrews_test(y3, "y3","both")
## 
## ################################ 
## # Zivot-Andrews Unit Root Test # 
## ################################ 
## 
## 
## Call:
## lm(formula = testmat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.8526 -0.7051 -0.0395  0.6218 24.0620 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.5411247  0.1711959  -3.161  0.00162 ** 
## y.l1         0.9700086  0.0064287 150.888  < 2e-16 ***
## trend        0.0068734  0.0014215   4.835 1.54e-06 ***
## du          -0.8168274  0.1928609  -4.235 2.49e-05 ***
## dt          -0.0019936  0.0006924  -2.879  0.00407 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.258 on 994 degrees of freedom
##   (1 observation effacée parce que manquante)
## Multiple R-squared:  0.9992, Adjusted R-squared:  0.9992 
## F-statistic: 3.135e+05 on 4 and 994 DF,  p-value: < 2.2e-16
## 
## 
## Teststatistic: -4.6653 
## Critical values: 0.01= -5.57 0.05= -5.08 0.1= -4.82 
## 
## Potential break point at position: 478

summary_y4 <- zivot_andrews_test(y4, "y4","both")
## 
## ################################ 
## # Zivot-Andrews Unit Root Test # 
## ################################ 
## 
## 
## Call:
## lm(formula = testmat)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -4.695 -0.837 -0.077  0.635 49.611 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.891355   0.305102  -2.921  0.00356 ** 
## y.l1         0.926948   0.010361  89.466  < 2e-16 ***
## trend        0.003445   0.001683   2.047  0.04090 *  
## du           3.228448   0.524691   6.153 1.10e-09 ***
## dt           0.009953   0.002304   4.320 1.72e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.42 on 994 degrees of freedom
##   (1 observation effacée parce que manquante)
## Multiple R-squared:  0.9984, Adjusted R-squared:  0.9984 
## F-statistic: 1.595e+05 on 4 and 994 DF,  p-value: < 2.2e-16
## 
## 
## Teststatistic: -7.0508 
## Critical values: 0.01= -5.57 0.05= -5.08 0.1= -4.82 
## 
## Potential break point at position: 299

Nous observons pour le test 1 une statistique de test supérieure aux valeurs critiques, on ne peut donc pas rejeter Ho et dire que la série est stationnaire avec une rupture. Cela conserve le fait que la série est une marche aléatoire pure.

Pour le test 2, on observe une statistique de test légèrement inférieure à la valeur critique en 5% en général donc on peut considérer potentiellement une rupture en 499.

Pour la série 3, on observe une potentielle rupture en 299. La série semble non stationnaire sauf s’il existe au point 299, rupture statistique significative. La statistique de test est inférieure à la valeur critique à 5% en général.

Pour la série 4, on observe potentiellement deux ruptures, une en 299 et une en 499, ce qui est également observé par la statistique du test inférieure aux valeurs critiques à 5% en général.

Les valeurs varient selon la marche aléatoire générer.

Le résumé ainsi obtenu des quatres séries est:

##   Series Test.Statistique Valeur.critique.1. Valeur.critique.5.
## 1     y1        -4.827725              -5.57              -5.08
## 2     y2        -6.139205              -5.57              -5.08
## 3     y3        -4.665259              -5.57              -5.08
## 4     y4        -7.050771              -5.57              -5.08
##   Valeur.critique.10. Point.de.rupture
## 1               -4.82              503
## 2               -4.82              498
## 3               -4.82              478
## 4               -4.82              299

4) Test pour Johnson & Johnson

Nous allons maintenant réaliser le test de zivet et andrews sur la data filtré de johnson et johnson.

Il peut être intéressant d’utiliser ce test afin de déterminer des potentielles ruptures dans une série.

## 
## ################################ 
## # Zivot-Andrews Unit Root Test # 
## ################################ 
## 
## 
## Call:
## lm(formula = testmat)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -12.4679  -1.2077   0.0247   1.4128  11.2114 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.966202   0.958805   4.137 4.57e-05 ***
## y.l1         0.797946   0.044296  18.014  < 2e-16 ***
## trend        0.030550   0.009281   3.291  0.00111 ** 
## y.dl1       -0.126280   0.063932  -1.975  0.04915 *  
## y.dl2       -0.003575   0.064165  -0.056  0.95560    
## y.dl3       -0.022461   0.062252  -0.361  0.71849    
## y.dl4        0.020673   0.061942   0.334  0.73880    
## y.dl5        0.091419   0.061495   1.487  0.13816    
## y.dl6       -0.087131   0.061300  -1.421  0.15624    
## y.dl7        0.191349   0.060417   3.167  0.00170 ** 
## y.dl8        0.178664   0.058768   3.040  0.00257 ** 
## du          -2.207872   1.008676  -2.189  0.02937 *  
## dt           0.139651   0.031366   4.452 1.20e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.171 on 302 degrees of freedom
##   (9 observations effacées parce que manquantes)
## Multiple R-squared:  0.9953, Adjusted R-squared:  0.9951 
## F-statistic:  5367 on 12 and 302 DF,  p-value: < 2.2e-16
## 
## 
## Teststatistic: -4.5614 
## Critical values: 0.01= -5.57 0.05= -5.08 0.1= -4.82 
## 
## Potential break point at position: 160

Nous avons réaliser le test sur la série représentant le prix de l’action et observons que aucune rupture significative ne semble se dégager. La statistique du test est supérieure aux valeurs critiques, on ne peut donc pas rejeter Ho.

## 
## ################################ 
## # Zivot-Andrews Unit Root Test # 
## ################################ 
## 
## 
## Call:
## lm(formula = testmat)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.4674  -1.5571   0.2017   1.6087  13.6416 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.004203   0.435534  -0.010 0.992308    
## y.l1        -1.034374   0.278140  -3.719 0.000239 ***
## trend        0.005830   0.002732   2.134 0.033683 *  
## y.dl1        0.894943   0.260989   3.429 0.000690 ***
## y.dl2        0.585029   0.240158   2.436 0.015430 *  
## y.dl3        0.374016   0.211199   1.771 0.077586 .  
## y.dl4        0.208421   0.183611   1.135 0.257228    
## y.dl5        0.143200   0.153894   0.931 0.352852    
## y.dl6       -0.090918   0.122558  -0.742 0.458767    
## y.dl7       -0.002327   0.088523  -0.026 0.979047    
## y.dl8        0.041763   0.059479   0.702 0.483131    
## du           3.856284   1.287869   2.994 0.002979 ** 
## dt          -0.220981   0.056876  -3.885 0.000126 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.438 on 301 degrees of freedom
##   (9 observations effacées parce que manquantes)
## Multiple R-squared:  0.1954, Adjusted R-squared:  0.1634 
## F-statistic: 6.093 on 12 and 301 DF,  p-value: 1.428e-09
## 
## 
## Teststatistic: -7.3142 
## Critical values: 0.01= -5.57 0.05= -5.08 0.1= -4.82 
## 
## Potential break point at position: 285

Dans un second temps, nous réalisons le test de Zivot et Andrews sur la série différencier des prix de l’action, qui est une série stationnaire.

Nous observons une statistique du test largement inférieure aux valeurs critiques, on peut donc rejeter Ho et dire que la série est stationnaire avec une rupture en 285.

Exercice 6 : Modeling the business cycle

Sur la base de toutes vos connaissances, proposez la spécification la plus appropriée pour modéliser la dynamique mensuelle de l’écart de crédit aux États-Unis. Les données sont disponibles sur le site web de la Fed de Saint Louis. Vous devez charger les taux d’intérêt annuels des obligations d’entreprise Moody’s Seasoned Aaa Corporate Bond et Moody’s Seasoned Baa Corporate Bond. L’écart de crédit est la différence entre l’indice Baa et l’indice Aaa.

Réalisez une analyse complète (saisonnalité, stationnarité, modélisation, vérification des résidus et prévision de l’écart de crédit). L’horizon de prévision est fixé à trois mois. L’étude devrait inclure des graphiques illustratifs, des commentaires détaillés et motivés, ainsi que des résultats pertinents.

Dans un premier temps, nous allons charger les données à partir de 1950.

getSymbols(c("AAA", "BAA"), from="1950-01-01",src = "FRED")
## [1] "AAA" "BAA"
AAA_monthly <- to.monthly(AAA, indexAt = "lastof", OHLC = FALSE)
AAA_ts=ts(AAA_monthly, start = c(1950, 1), frequency = 12)

BAA_monthly <- to.monthly(BAA, indexAt = "lastof", OHLC = FALSE)
BAA_ts=ts(BAA_monthly, start = c(1950, 1), frequency = 12)

spread_ts=BAA_ts-AAA_ts

On peut observer que le spread est une série ayant plusieurs tendances avec potentiellement une saisonnalité mensuelle. De plus l’ACF de la série semble décroître lentement, cela fait penser à un processus avec mémoire. De même avec la PACF, qui a un pic très fort au lag 0 ce qui explique la décroissance lente de l’ACF, toute la mémoire dépend du premier paramètre autorégressif.

On peut donc penser que cette série est non stationnaire.

Afin de se rendre compte de la saisonnalité et des différentes tendances, nous allons utiliser la fonction decompose.

Afin de pouvoir enlever cette saisonnalité, nous allons utiliser un filtre de moyenne mobile précedemment décrit.

On remarque sur le graphique que l’on a bien réussi à enlever la saisonnalité.

L’ACF décroit toujours lentement et le premier coef de la PACF n’a pas changé. On peut donc penser que la série n’est toujours pas stationnaire. On va vérifier cela avec un test de stationnarité, ici le test de KPSS.

## 
## ####################### 
## # KPSS Unit Root Test # 
## ####################### 
## 
## Test is of type: tau with 6 lags. 
## 
## Value of test-statistic is: 0.7785 
## 
## Critical value for a significance level of: 
##                 10pct  5pct 2.5pct  1pct
## critical values 0.119 0.146  0.176 0.216

On remarque ici que la valeur de la statistique est supérieure aux valeurs critiques, on peut donc rejeter Ho et ainsi dire que la série est non stationnaire.

Afin de rendre la série stationnaire, nous allons faire une différentiation avec la fonction diff afin de supprimer les tendances et ensuite de nouveau tester la stationnarité avec le test de kpss et également ADF.

spread_diff=diff(spread_ma,differences=1)
kpss_test <- ur.kpss(spread_diff,type='tau')
summary(kpss_test)
## 
## ####################### 
## # KPSS Unit Root Test # 
## ####################### 
## 
## Test is of type: tau with 6 lags. 
## 
## Value of test-statistic is: 0.0217 
## 
## Critical value for a significance level of: 
##                 10pct  5pct 2.5pct  1pct
## critical values 0.119 0.146  0.176 0.216
test=ur.df(spread_diff,type="none",selectlags="BIC") 
summary(test)
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression none 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.083091 -0.004704  0.000172  0.004873  0.083516 
## 
## Coefficients:
##            Estimate Std. Error t value Pr(>|t|)    
## z.lag.1    -0.10266    0.01195  -8.594   <2e-16 ***
## z.diff.lag  0.39393    0.03113  12.656   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01255 on 872 degrees of freedom
## Multiple R-squared:  0.1863, Adjusted R-squared:  0.1844 
## F-statistic: 99.81 on 2 and 872 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -8.5937 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau1 -2.58 -1.95 -1.62

Pour le test de KPSS, on observe une statistique de 0.02, inférieure aux valeurs critiques, on peut donc conserver Ho.

Pour le test d’ADF, on observe une statistique largement inférieure aux valeurs critiques. On peut donc rejeter Ho et dire que la série est stationnaire.

Le degré d’intégration de la série de spread est donc de 1.

Nous allons maintenant pouvoir estimer un model ARMA(p,q) pour notre série.

On observe que l’ACF se rapproche de 0 vers le lag 3 ou 4 et que la PACF se rapproche de 0 vers le lag 3. On va donc vérifier cette hypothèse en calculant les critères d’information pour des valeurs de p et q entre 1 et 4.

##           [,1]      [,2]      [,3]      [,4]
## [1,] -5131.692 -5129.987 -5130.191 -5133.521
## [2,] -5133.902 -5188.773 -5196.916 -5225.748
## [3,] -5137.494 -5212.798 -5206.045 -5222.845
## [4,] -5144.154 -5206.067 -5204.906 -5204.564
##           [,1]      [,2]      [,3]      [,4]
## [1,] -5117.328 -5110.836 -5106.252 -5104.794
## [2,] -5114.751 -5164.834 -5168.189 -5192.233
## [3,] -5113.555 -5184.071 -5172.530 -5184.542
## [4,] -5115.427 -5172.552 -5166.604 -5161.473
##           [,1]      [,2]      [,3]      [,4]
## [1,] -5126.201 -5122.666 -5121.040 -5122.539
## [2,] -5126.581 -5179.622 -5185.935 -5212.936
## [3,] -5128.343 -5201.816 -5193.233 -5208.203
## [4,] -5133.172 -5193.255 -5190.264 -5188.091
## [1] "Optimal AIC - p: 2 q: 4"
## [1] "Optimal BIC - p: 2 q: 4"
## [1] "Optimal QIC - p: 2 q: 4"

On va donc modéliser notre série sous la forme d’un ARMA(2,4).

## 
## Call:
## arima(x = spread_ts, order = c(2, 1, 4))
## 
## Coefficients:
##           ar1      ar2     ma1     ma2     ma3      ma4
##       -0.2931  -0.5656  0.6506  0.6362  0.1151  -0.1252
## s.e.   0.2757   0.1127  0.2729  0.1509  0.0594   0.0351
## 
## sigma^2 estimated as 0.009606:  log likelihood = 800.59,  aic = -1587.18
## 
## Training set error measures:
##                        ME      RMSE        MAE        MPE     MAPE     MASE
## Training set 0.0003114853 0.0979554 0.05997818 -0.1977236 5.759479 0.971493
##                     ACF1
## Training set -0.00193067

On observe ici un RMSE plutôt faible, cela montre que le model a bien approximé notre série.

Nous allons maintenant procéder aux tests de qualité: ditribution normale des résidus, absence d’autocorrélatio, homoscédasticité et stationnarité.

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals_model
## W = 0.83844, p-value < 2.2e-16

On observe que notre histogramme semble correspondre à celui d’une loi normale mais les quantiles extrèmes de nos résidus s’éloignent de ceux d’une loi normale. De même, le test de shapiro affiche une p-value très faible donc on rejette le fait que les résidus se comportent comme une loi normale.

## 
##  Box-Ljung test
## 
## data:  residuals_model
## X-squared = 17.46, df = 10, p-value = 0.06479

La p-value du test de Ljung Box est supérieure à 5%, on peut donc conserver le fait qu’il n’y ait pas d’autocorrélation dans les données. Ce fait peut également être vérifier avec l’ACF qui converge rapidement vers 0.

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression none 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.48996 -0.03794 -0.00332  0.03596  0.90540 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## z.lag.1    -1.004253   0.047638 -21.081   <2e-16 ***
## z.diff.lag  0.002329   0.033653   0.069    0.945    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.09818 on 883 degrees of freedom
## Multiple R-squared:  0.501,  Adjusted R-squared:  0.4998 
## F-statistic: 443.2 on 2 and 883 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -21.0809 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau1 -2.58 -1.95 -1.62

En réalisant le test d’ADF sur les résidus, on trouve une statistique de -21, largement inférieure aux valeurs critiques. On peut donc rejeter Ho et dire que la série des résidus est stationnaire.

Le dernier test à réaliser est celui de l’homoscédasticité des résidus.

On observe que la PACF des résidus au ² converge rapidement vers 0 après le lag 1. On va donc estimer un model AR(1) pour le carré des résidus.

## 
## Call:
## ar(x = residuals_2, order.max = 1)
## 
## Coefficients:
##      1  
## 0.1571  
## 
## Order selected 1  sigma^2 estimated as  0.001409

Le coef de notre AR est proche de 0, on peut donc dire que les résidus sont homoscédastiques.

De plus nos résidus semblent être bien distribués et à variance constante.

Pour conclure, on peut dire que nos tests de qualité permettent d’assurer la validité du model.

On va finir par la prédiction de valeurs sur 3 mois.

##                80%       95%
## Dec 2023 0.8800539 0.8135623
## Jan 2024 0.7945037 0.6823941
## Feb 2024 0.7379131 0.5953966
##               80%      95%
## Dec 2023 1.131265 1.197757
## Jan 2024 1.218064 1.330173
## Feb 2024 1.276353 1.418870

Nous avons afficher les tableaux des bornes des intervals de confiance